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I. INTRODUCTION 

Due to its unique band structure, the charge carri- 
ers in graphene are massless Dirac fermions which can 
cross high potential barriers with ideal unity transmis- 
sion coefficient (the Klein paradox^. This ensures a 
very effective escape channel from a trapping potential 
thus making it hard for conventional Dirac electrons to 
be localized within graphene based QDs. Within a finite 
spatial region defined by sharp potential profile*"^. To 
overcome this difficulty, and trap the electrons for suffi- 
ciently long time, we propose to dress the electrons with 
circularly polarized photons, thus providing them with 
an effective mass^l. The localization is demonstrated 
in a cylindrical quantum dot (QD) formed in monolayer 
and bilayer graphene by antisymmetric potential kink. 
Conventionally the measure of localization are charac- 
teristic resonances in the electronic density of states^'. 
The dynamical gap is studied semi-classically using Flo- 
quet's theorenP23 We present a fully quantum mechani- 
cal model, which is based on dressing electrons in mono- 
layer and bilayer configurations. Our calculations show 
that the dressing not only opens up a dynamical gap in 
the energy dispersion but also renormalizes the Fermi ve- 
locity and inter layer coupling coefficients. In the bilayer 
configuration, the dressing tunes the gap. That is, it 
can either close or open the gap, depending on the po- 
larity of the potential kink and the direction/degree of 
the polarization. The resulting confined electronic states 
should have similarities with the surface states of topo- 
logical insulators. Their energies are located inside the 
energy gap and the wave functions decay away from the 
interface of the kink potential. These topological states, 
with the carriers propagating along the potential kink, 
are expected to be robust with respect to the effects 
of disorder^!. The fully localized states are mixed with 
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the quasi-bound states above the energy gap which can 
effectively carry away the energy. Conventional linear 
response spectra (proportional to the density of states) 
provides limited information about them due to the large 
broadening caused by their short lifetime^ . We propose 
to utilize femtosecond nonlinear spectroscopy in order to 
study their dynamics. We shall use a four-wave mixing 
technique known as photon-echo^ The mixed signal is 
heterodyne detected in the direction of — ki + k2 + k3 as 
shown in Fig. [3] The photon echo is known to be able 
to eliminate the inhomogeneous broadening due to im- 
purities, and allow us to focus on the intrinsic lifetimes 
of the electronic states. We further discuss signatures 
of the dynamic gap on the two-dimensional (2D) spec- 
tra. There is yet another peculiar characteristic of local- 
ized Dirac electrons. As in metals, they are dynamically 
screened, leading to small Coulomb interaction between 
them. For small QD this leaves Pauli blocking as the pri- 
mary source of the nonlinear signaP. This allows us to 
calculate it as a sum-over-states (supermolecule) formal- 
ism. We can further simplify the signal interpretation by 
switching to the quasiparticle picture. Those are given 
as deviation from ideal bosons^ for which the nonlinear 
signals vanishes. We are able to consider only excited 
states absorption Liouville pathways without contribu- 
tion from the ground state bleaching and excited states 
emission. This interference reduction provides relatively 
simple interpretation of the 2D spectra. The short lived 
states can be visualized via the coherences (off-diagonal 
cross resonances) with those fully localized. We employ 
visible light to map the QD interband transitions onto 
the 2D spectra. Finally we briefly discuss the effect of 
the Coulomb induced exciton scattering based on nonlin- 
ear exciton equations!^. Possible formation of biexciton 
molecules is demonstrated. 

The outline of this paper is as follows. In Sec. [TTJ we 
present the model Hamiltonian for graphene irradiated 
with a circularly polarized electromagnetic field. Section 
|III| is devoted to dressing of electrons in bilayer graphene 
and a derivation of the eigenstates. We deal with the 
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trapping of the dressed electrons within a QD in Sec. 
IV Section [V] presents the absorption and correlation 
spectrum for dressed electrons in a QD. We present nu- 
merical results in Sec. I VII and conclude in Sec. I VII with 
a summary of our results. 



II. DRESSED ELECTRONS IN FREE STANDING 
GRAPHENE. 

The electronic Hamiltonian of graphene irradiated 
with an electromagnetic field may be expressed as^: 



H — Ho + "Hi + H2 



H = hioo%a Q 



(J+a + (T-aJ 



Ho 



Hi = hvpcr ■ k + IV(x, y) 



h^ l a\a l + —j= ({(7+ + a-)(ai + at)) 



(1) 

(2) 
(3) 
(4) 



Here, Ho describes the Jaynes-Cummings modeP^l with 
ao being the annihilation operator of a single mode cir- 
cularly polarized optical field with frequency uiq and ./Vo 
photons in the mode. Each of them carries the energy 
hajo- cr± — (cr x ± icr y )/2 are raising and lowering oper- 
ators for z— projection of the electrons pseudo-spin. In 
matrix representation these are Pauli matrices. Wo/\/No 
is the electron-photon coupling, a quantum mechanical 
analogue of the classical rotational motion caused by the 
circularly polarized wave. 

Hi describes conventional Dirac HamiltoniarpJ of 
graphene with Fermi velocity vf = 10 6 m/s; k is the 
wave vector measured from one of the K points, V(x,y) 
is an external QD confining potential; X is the identity 
matrix. H2 describes the rest of the optical modes lat- 
ter used to probe the dressed states by four-wave mixing 
process. 

The Hamiltonian Hq may be diagonalized in a straight- 
forward way^ in the following basis: 



_ ( \^+,N ) 



|V>±,aO = cos0|±,iVo) ±sin#F,iVo±l) 



in 



COS ( 



sm = 



A',, 



f2jv — fiwo 

2/ 



n No = nuo + w^N + iyN 



(5) 
(6) 

(7) 
(8) 



Here, the direct product state |±,iV ) define the un- 
coupled electron with pseudo spin up (+) or down (-) 
and the optical mode with iVophotons. Eq. ^ defines 




FIG. 1. Bilayer graphene structure and renormalized cou- 
pling coefficients. 



the dressed electron states. In the basis of Eq. ([5]), the 
Jaynes-Cummings Hamiltonian assumes the form 



(lpN \H \^No) 

(lp+,N \Ho\^+,N ) (lp+,N \Ho\lp-, No ) 
(l()- lNo \Ho\ll>+,No) {^-,N \H \^-,N ) 



(9) 



N hujQ + E g /2 







N htj - E g /2 

= INoHloo + (E g /2)a 3 . 

The first term is a constant, and may be omitted. 

The remaining Hamiltonian matrix elements are cal- 
culated in Appendix [B] yielding 

(ip No \ni\Tp No }=M F <T-k+TV(x,y) (10) 

(^N \H 2 \4>N ) = (11) 

m 



00 VV 



where vp = vf cos 2 <f> is the renormalized Fermi velocity 
and Wi = Wi cos 2 cj> are the renormalized couplings to 
the probing optical modes. 

In the absence of a potential (V(x,y) — 0), the eigen- 



values of Ho + Hi are ±J (hV F k) 2 + (E g /2) 2 and the 
corresponding eigenfunctions are: 



#+(fc) = e lkr 
*+(fc) = e 4kr 



cos (ttfe/2) 
e t/3k sin(a fe /2) 

sin (a k /2) 
-e^" cos (ak/2) 

Tan/3 fc = k y /k x ; 

Tanafe = 2hvpk/ E g . 



(12) 

(13) 
(14) 



III. DRESSED ELECTRONS IN BI LAYER GRAPHENE 

Starting with Eq. (38) of the review article of Castro 
Neto, et al.^, and applying the procedure of Appendix 
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\K\ the Hamiltonian which describes the dressing of the 
electrons in the bilayer (Bernal stack) can be expressed 
as 



H = 



Hu H 
H21 H2 



12 



(15) 



Here Hu = H^Yi) an d #22 — UiVz) describe the elec- 
trons on the first and second graphene layers respectively. 
Those layers may experience different potential profiles 
(Vi.,2) entering Eq. (J3J). The interlayer coupling is de- 
scribed by the off-diagonal block matrices as: 

H12 = 7iO"_ + 3j 3 a(k x - ik y )a + (16) 
W§ W i < \ t 1a 



'N, 



H 2 i = li<J+ + 3j 3 a(k x + ik y )<T_ (17) 
("era. + a- 



= V-4 + -7#( (7 + + <r-)(at + 4) 



i=i 



Here, we have introduced effective electron-ph oton cou- 
pling matrix elements W^/y/Nl = —3ejj a y // 2n/uji £lh, 
where a = 1.42 A is the carbon-carbon distance within 
a layer. Additionally, we have 70 = 2.8 , eV is the 
nearest-neighbor hopping energy within the layer (A\ ^ 
B\, A2 ^ B2). Fermi velocity can be expressed in terms 
of the above parameters as hvp = 37oet/2. 71 = 0.4 eV is 
the inter-layer hopping energy between atoms of type A: 
(A\ ^ A2). 73 = 0.3 eV is the inter-layer hopping energy 
between atoms of type B: (Bi ^ B2). 74 = 0.04 eV is the 



inter-layer hopping energy between atoms of type B and 
A: (Ai ^ B2,Bi A 2 ). Electronic couplings between 
various atoms in the bilayer graphene are shown in Fig. 
[T] The double layer can be regarded as corresponding to 
as 71 = 73 = 74 = 0. 



In the dressed state basis of Eq. (C2), the diagonal 



blocks are given by the results of the previous section as 

(lpN \Hll,22\i>N ) (18) 

= hv F a ■ k + (E g /2)a 3 + IVi,»(x, y) 
+!N hu} + )hw i a\a i + — =(0+ + 0-)( a i + a}) . 



The off-diagonal blocks may be derived from Eqs. (16 1, 
(17) and Appendix [C] to become: 



(ipNolHul^No) = 7i<7- + 3j3a(k x - ik y )(T+ 

+74^3 + ~JW^+ + <J -^°' 1 + °<) 



i=l 



(19) 



(20) 



(■>Pn \H2i\iI'n ) = 71°"+ + 37 3 a(fc K + ik y )a- 

+74^3 + J2 V 7k ( - (J + + a -^ a * + a ^ ' 

where the renormalized model parameters are 712 — 
71,2 cos 2 <j), 74 = (W / 3 /2) sin 20 and Wf = Wf cos 2 4>. For 
the purpose of further discussion, it is convenient to local- 
ize Hq+Hi as we did in a preceding section for monolayer 
graphene. The corresponding matrix elements are 



Single layer: 



Bilayer: 



Vi{x, y) + (E g /2) l^ a(k x + ik y 



3 z. 



j a(k x - iky) V\ (x, y) - (E g /2) 



{ip No \{B2A2A 1 Bi\H + H 1 \B 1 AiA 2 B2)\ilJ No 
( Vi(x,y) + (E g /2) ^%a(k x + ik y ) 74 Slza(k x - ik y ) ^ 



3 x 



loa(k x - iky) Vi(x, y) - (E g /2) 



7i 



"74 



74 



\ 37 3 a(fc x + iky) 



71 
-74 



V 2 (x, y) + (E g /2) §7bo(ftx - iky) 



3- 
2 



j a(k x + iky) V 2 (x, y) - (E g /2) ) 



(21) 



(22) 



This implies that the dressing of the Dirac electrons in 
bilayer gives 

• renormalized interlayer coupling coefficients, which 
are denoted by tilde., 

• broken the symmetry between the sub-lattices 
(Ai, Bi; A 2 , B 2 ) of each of the layers. Measure of 
the broken symmetry is (E g /2), 



• broken symmetry between the sub-lattices 
(Ai, B2; A2, Bi) belonging to different layers. A 
measure of the broken symmetry is 74. 

The corresponding eigenvalues for constant potentials 
(V\, V2) are shown in Fig. [2] for chosen values of the pa- 
rameters. We first focus on the largest interlayer coupling 
7x and neglect the rest of the coupling (Figs. [5Ja)). The 
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four bands are given by 



{2E + V 2 + V l ) 



(23) 



E 2 g + (Vi - V 2 f + 9a 2 fc 2 7 2 + 2 7l 2 ± 2^ (Vi - V 2 ) 2 + 9a^ 2 7 2 ) + 7? (9a 2 fc 2 % 2 - 2E g (Vi 



v2)) + ii ■ 



On its own, E g opens a gap in the bilayer spectrum sim- 
ilar to the monolayer (Fig. [2ja.l)). The gap may be 
opened by applying a potential difference between the 
layers (Vi ^ V 2 in Fig. (2ja.2)). The combined effect 
of the potential difference and E g > can either widen 
V2 — Vi < or shrink V% — V\ > the gap compared 
with the gap induced by the potential difference itself 
(Fig. [2];a.3)). We observe that when 2E g =V 2 -Vi, the 
gap closes (Fig. [2]ja.4)). Inclusion of the rest of the cou- 
pling breaks the symmetry between k x and k y , as follows 
from Fig. [T] The analytical form of the energy bands, 
although possible, is too large to be presented here. The 
energy bands are shown in Fig. [2jb-d). 



IV. DRESSED ELECTRONS CONFINED IN A QD 

Let us now turn to the problem of trapping dressed 
electrons within a QD. Since the confining potential 
Vi(x, y) is radial it is convenient to work with cylindrical 
coordinates x — r cos 9 y — r sin 8. This ammounts to 
the following substitutions: 



k x id x , ky 



d x ± id y = e ±ie ( d r ± "-d t 



(24) 



Thanks to the potential radial symmetry the Hamilto- 
nians in Eq. (21) and (22) commutes with the angular 
momentum operator L z = I{xk y — yk x ). Here we have 
neglected symmetry breaking contributions to the Hamil- 
tonian (73,4 <C 71). Therefore, we may seek solution of 
the Dirac equation in the form: 



l*m(r)> = 



(25) 



Here, the projection of the angular momentum has eigen- 
values to = ±1/2, ±3/2, ±5/2,.... Substituting Eqs. 



(25) and (j24h into (22), we obtain the following system 



of ordinary differential equations 



( V 1 (r) + (E g /2) f 7oa (a r -^) 
-§7oa(a r , + ^ 2 ) V l {r)-{E g /2) 








71 







7i 

V 2 (r) + (E g /2) 

3 ~ / o m— 1/2 

§7 a(d r 



3~ m+1/2 



) 

) V 2 (r) - (E g /2) J 



^m,B x { r ) 

^ m , Al (r) 



E 



^ m .B l {r) 

i> m ,A 2 {r) 
ip m ,B 2 (r) 

(26) 



First, let us consider the case when there is no coupling 
between the graphene layers. Assuming that the solution 
of Eq. ( 26 ) has the form of y/ripm (r) in the regions of 



constant potential, we obtain 



ip m ,B 1 (r < ) 
Tp m ,Ai( r <) 



-0m,-Bi( r >) 
, 4' m ,A 1 { r >) 



( -Jm 1/2 [ ;;-,,„ 



2r<V^-(g s /2)^ \ 



. 2r <y /E2-(E g /2) 
J \m+l/2\ \ 3^ 



( ff,W 



(27) 



= B 



H 



|m+l/2| 



2r > i/(E-V 1 y-(E g /2y } \ 
/ 



\tn— 1/2 1 I 370a 
<l) ( 2r >y /(E-V 1 )*-(E g /2)* 



370 a 



(28) 



The Bessel function form of the wave function inside of 
the QD (Vi = 0, r = r< < R) is dictated by the fact 
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FIG. 2. Dotted/solid curves represent electron dispersion of double/bilayer graphene. Panels (a) assume that 71/70 = 0.4 
and the remaining interlayer coupling parameters are zero. Panels (b) also introduce the effect of 73/70 = 0.3. Panels (c) 
demonstrate the effect of 74/70 = 0.3. Panels (d) illustrate the combined effect of 73 and 74. In columns (1) we have dressing 
induced i? 9 /7o = 1/2 with no potential difference between the layers. In columns (2) we potential difference between the layers 
(Vi — V"2)/7o — 1/2 with no electron dressing Eg/yo = 0. Columns (3,4) illustrate their combine effect; (3) correspond to 
Eg/^o = 1/2 and (4) has E g /-yo — 1. Blue (green) curves show a section of the energy along k x (k y ) directions. 



that the wave function must stay finite at r = 0. Outside 
the QD (Vi > 0, r — r> > R) the wave function must 
describe the outgoing wave at large distances (r> 3> R). 
We, therefore, took it to be the Hankel function of first 
kind. At the boundary of the dot the wave function must 
be continuous. The energies E m ^ n of the quasi-stationary 
states inside of the QD are obtained by solving the fol- 
lowing equation: 



H, 



\rn— 1/21 ( k > R ) _ J\ m— 1/2| 

(k<R) 



H, 



(i) 



'|m+l/2| "W^l 

where we have introduced following notation: 

_ 2Ry/(E V,r (E g /2f 

K > — 5^ 

370 a 

_ 2Ry/E* - {E g /2Y 

K < — 5^ 

570 a 

Those are shown in Fig. [3] for several chosen values of 
Eg. The long living solutions Im[E] ss can be obtained 



(29) 



analytically by noticing the following identities for the 
Hankel function in the limit z«l, 



H^(z) = J n (z)+iY n {z) 



(30) 



1 



r(n + i) 

It is clear from the above equation that when E = V\ ± 



Eg/2 the left hand side of Eq. (29) vanishes. Therefore 



the real energies of the QD correspond to the zeros of the 
Bessel function with 



J\m—1 



/2| 



2R+/E n 



= 



(31) 



The splitting of the central peak in the density of states 
(DOS) by the electron dressing should be readily acces- 
sible in optical experiments. This will be the subject of 
the following section. 
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FIG. 3. Energy levels (Blue dots) of a QD formed in a 
graphene monolayer. Green dots correspond to dressed Dirac 
electrons. The confining potential is Vi = 20, and the gap 
E g = f0 (in units of 2R/3y a). 



Now let us briefly comment on the QD in the bi- 
layer graphene. Since the dressing of the electrons in 
the bilayer allows wide manipulation of the gap (See 
Fig. [2]ja.3)) one can dynamically form conventional QDs. 
That is those which support infinitely living electronic 
states. There are two possible schematics shown in Fig. 
In one case the dressing and the substrate induced gap 
work concurrently and the whole realization of the QD 
is almost identical to the one we had discussed for the 
single layer. The only difference is the non-homogeneous 
potential which forms the dot is the potential between 
the graphene layers. Such QD can be readily realized by 
screening of the substrate potential inside of the QD. In 
the other schematic, the combination of the dressing and 
the substrate potential closes the gap inside of the QD. 
The electrons become trapped by the gap outside of the 
QD. In both cases the analytical solution of Eq. (26) 



may be obtained bv following the procedure outlined in 
RefPH 



V. ABSORPTION AND PHOTON-ECHO SPECTRA OF 
DRESSED DIRAC ELECTRONS IN A SINGLE QD 

In this section we investigate several linear and non- 
linear optical techniques which allow to probe the details 
of the electronic structure calculated above. The linear 
absorption mostly reveals the long living states (lm[i?] <C 
Re[_B]) which show up a narrow res onance and directly 
reflects the structure of the DOS^l Nonlinear photon- 
echo^ signal ( x^ 3 H~ ki + k2 + k.3)) will be designed to 
reveal the other short living states. The schematic of 
the heterodyne detected four wave mixing experiment is 
shown in Fig. [4] By probing the coherence between the 
electronic states in the QD, the technique can reveal the 
energy of the short living electronic quasi-bound states. 
We shall restrict the discussion to singly excited states, 




FIG. 4. Schematic illustration of the photon-echo technique 
designed to study the exciton scattering dynamics in graphene 
based QDs. 



thereby neglecting underlying many-body effects. This 
allows for a conceptually simpl e description in terms of 
the many body eigenstate d 11 * 16 *. 

It is also possible to study several excitonic states of 
the QD and their dynamics by applying special nonlin- 
ear measurements. We shall also make use of the double 
quantum coherence^ (x^C^-i +^-2 ~ k^)) technique in 
order to observe many body effects in biexciton manifold 
of graphene QD. Con ceptu ally the approach is similar to 
that described in Ref! 13 * 17 ( However the graphene based 
QD has large screening of the Coulomb interaction be- 
tween excitons. Thus it can be safely omitted in the 
following discussion. Recently the collinear version of 
X^ 3 ) technique based on phase-cycling gained popularity 
in QD studies of terahertz regime 12 . Its application to 
graphene will be reported elsewhere. 

Let us first define the effective single particle diagonal- 
ized Hamiltonian in the QEP^ 



ml,nl ^ m iC n i 

ml,nl 

+ E 



(32) 



m2,n2 



whose matrix elements are obtained from Eq. ( 29 ). Each 



subscript is a composite of two indices describing angular 
(m) and radial (n) quantum numbers ml(2) = {m,n}. 
Here we have partitioned the electronic states into oc- 
cupied (n < 0) and unoccupied (n > 0) in the ground 
state obtained by setting up the chemical potential to 
fi = NqHluq. The electrons in the unoccupied state can 
be created by action of & ml = |* m , n >o) ® (^m,n<o\ 
operator on the ground state. Its hermitian conjugate 
c m \ removes the electron from that state. Similarly 
the second term of Eq. ( 32 ) describes the creation 

and annihilation d m i of 



/t 



h m1 — \*m,n<0) ® (^m,n>0 

the electrons in the originally occupied states (holes). 
The second term in Eq. (32 1 is just hermitian conju- 
gate of the first. In the notation above symbol <g) stands 
for element-wise multiplication of the vectors. Since the 
dressing circularly polarized CW mode has been already 
incorporated in Eq. ( 32 1 we only have to explicitly 
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treat the interaction with the narrow pulsed time or- 
dered incoming and detected modes. This is given by 
the effective^ interaction Hamiltonian in rotating wave 
approximation: 



Hut(t)= ^™l,™2 rf m2C m l +H.C 

i,ml,m2 J=i 

(33) 

Here £*' stands for the left(right) ± polarized com- 
ponent of the incoming or detected mode electric field 
amplitude. The dipole moments of transitions (in units 
of Wi/(£ fVNl)) are: 



Direct diagonalization of the above Hamiltonian ( 38 ) in 



/Ci,m2 = J dr(*„ a (r)|re ±iy |v]/ m2 (r)) (34) 

Note that we are still in the single electron-hole repre- 
sentation, not yet in the many body exciton/hole repre- 
sentation. Therefore we do not need the envelop function 
to define the transition moments as in RefP^. 

The next step is to bring the Hamiltonian in Eq. 



(32 1,(33) into the excitonic form. Using the method first 
proposed in RefP^ we define the electron-hole pair anni- 
hilation operators (not to be confused with exciton oper- 
ators) as: 



fit, = Br, 



(35) 



where we used composite index of m — {ml, m2} Since 
we are interested in the third-order response the com- 
mutator of the above operators may be truncated at 
quadratic order: 



B n 



, B\ = 5 m , n - 2 22 5, 



S" t R 
. j, 1 'q 



(36) 



p.q 



where S m , n = <5 m i,ni£ m 2,n2- The tetradic matrix 8 m ,n- P ,q 
(phase-filling factor) is responsible for the deviation from 
the boson statistic of the pair operators, and steams from 
the fermionic nature of its constituents: 



25 



m,n;p,q 



(37) 

= < 5ml,gl < 5m2,p2<5nl,pl<5n2,g2 + ^ml,pl<5m2,g2^nl,(jl^n2,p2 

In the basis of electron-hole pairs the Hamiltonian in 



Eqs. (31), (32) becomes Frenkel-like if truncated up 
to forth order (valid for third order response with two 
excited electron- hole pairs): 

H = ^ E m B] n B m + - ^ (E m + E n ) B^B^BjnBn 

(38) 

n int (t) = ^^^-f 1 K4 + H.c (39) 

i,m j=± 



order to find the exciton/biexciton manifolds is difficult 
and non-equilibrium Green's functions for the single and 
double electron-hole pairs are used instead. If one ne- 
glects the nonlinearities caused by the Pauli exclusion 
those retarded Green's functions are defined as: 



Gel, 



G el (r) = -iO{r)e- iE ^ 



e2 



(40) 
(41) 



Here 8(t) is the Heaviside function and the time between 
two consecutive pulses is denoted as Tj = t i+ i — ij. Note 
that in order to be retarded the Green's functions must 
contain the energies with Im[£' e ] < 0. We also adopted 
the notation E e i = E m \ + E m 2- We shall also need their 
Fourier transforms with respect to the time delays: 



GelM = 



1 



LO - Eel 
1 



LO — E e i — E e 



(42) 
(43) 



In the above Green's functions the biexciton energies (the 
poles of Eq. ( |43| ) is simply a sum of the exciton energies. 
The nonlinear signal from such system vanishes since it 
represents a collection of harmonic oscillators. The effect 
of Pauli exclusion in Eq. ( 38 ) is usually incorporated by 



tetradic exciton scattering matrix: 



r 



e4,e3;e2,el 



(io) — 5 e l,e3;e2,el{G 1 )e4, 



<;3 



(44) 



Which carries all the information about underlying non- 
linearities. Coulomb interaction can be incorporated by 
solving Bethe-Saltpeter equation as in RefPHU]. 

The photon echo signal can be recast in terms of non- 
interacting Green's functions as well as the scattering 
matrix as: 

^f+i 3 2 f k3 (^3,r 2 =0, Wl ) (45) 

= 2Re Y Me3Me2 3 Vef Vet G* 3 (~Wi )G e4 (^ 3 ) 
el,e2,e3,e4 

x re4,e3;e2,el(w3 + E e 3)G e 2,el (<^3 + E e 3) 

Even though the above two expression formally gives 
the nonlinear signal, it is hard to analyze. However it is 
convenient for numerical simulations due to expandabil- 
ity of the scattering matrix into the domain where the 
coulomb interaction may play its role. The detailed form 
of the scattering matri x which involves the Coulomb in- 
teraction is given in Ref! 13 * 16 ^. An alternative approach 
to derive the signal by using double-sided time ordered 
Keldish diagrams shown in Fig. [5] The diagrammatic ap- 
proach (also known as " sum over states" ) can answer one 
of the fundamental question whether the nonzero scat- 
tering matrix is sufficient to calculate a nonlinear signal. 
The answer to that question is not trivial due to large 
number of interfering terms in Eqs. (45) and (Dl). The 



diagrams were constructed by blocking the consequent 
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FIG. 5. Feynman diagrams for the photon echo technique 
k s = — ki + k 2 + k 3 . 



VI. NUMERICAL RESULTS AND DISCUSSION 

The main advantage of the Pauli blocking description 
of excitons is its simplicity. For a model with N singly 
excited electronic states, we only need to consider N dou- 
ble excited states compared with N(N — 1) in the case of 
Coulomb scattering induced biexcitons. We note that the 
same number of doubly excited states N(N — 1) are al- 
lowed in a simple boson harmonic model. This allows for 
better tracking of pathways interference and resonances. 
In this section, we shall use it to classify the off-diagonal 
resonances in the 2D photon echo spectra in accordance 
with the short living excited states of the QD. This will 
be compared with the linear absorption spectrum which 
is proportional to the single excited electronic density of 
states and is given by the main diagonal of the 2D spec- 
tra. We shall demonstrate the improved resolution of 
those short living states via the coherent response with 
the long living excitation. First, we assume the a model 
of two single excited states (fi — P2,-E 2 ). This 

leads by the Pauli exclusion principle to the single dou- 
ble excited state (fii2, Ex+E^). The photon echo signal 
contains three distinct pathways: ground state bleaching 
(GSB Fig. gb)), excited state emission (ESE Fig. gd), 
and excited state absorption (ESA Fig. |5{a,c). Those 
are given by 



double excitation of the same electron-hole pair. The 
nonlinear signals can be ext racted from the diagrams by 
the rules stated in Ref! 11 * 16 ! In our case the photon echo 
signal is obtained via the diagrams in Fig. [5| 



= Re 



^ -ui-E. 



i'-ki+ka+k;, (^3,T2,-Wl} 

,.3'1 ,,*,& ,,*>J' 3 ,,J'4_ 
rm Pro r m ' Pro 



(46) 

-ir 2 (E m -E*J 



&3 — E„ 



EL 



E n 



*J2 * j3 7"4 -ir 2 {E m -E* ml ) 
Pra' rm r m ,e 



-uji — E„ 



^3 — E„ 



-Re 



m,m" 



-UJl + E; 



- E' m + E' m 

*,j2 *J3 ji 
rm Pm" Pro 

W3 — E rn 



+ 



*,i2 + 03 j'4 
Pro" M TO Jd Mro" e 



-iT 2 (E m -E^,) 



— -E„ 



This signal would vanish if it were not for the Pauli block- 
ing which prevents m to be equal to ml . 

In Appendix [D] we have analyzed an alternative form 
of the four-wave mixing known as double-quantum coher- 
ence. This signal vanishes identically despite the Pauli 
induced scattering since the exited state absorption path- 
ways are fully compensated by their ground state bleach- 
ing and exited state emission counterparts. Therefore the 
double-quantum-coherence can be readily used as a mea- 
sure of the Coulomb interaction strength, and screening. 



Sgsb( w 3,£2,wi) = Re 2^ 7— _ F * 



^(-^-EDius-Ej) 



SEffi(w3,T 2 ,a;i) = Re ^ 



(47) 



(48) 



= - Re E 



- Re E 



, J= i (""I - £*) (W3 - #i) 

S^lCwa.Tn.W!) (49) 
(-w x - E*) (w 3 - ^ - ^ + £*) 

2|,,.|2--n.(B*-B?) 



-wi - (w 3 - Ei - £y + E*) 



Clearly, when Pauli blocking is neglected, we have a col- 
lection of damped oscillators and at t 2 — the signal 
disappears. It would also vanish if we assume that there 
is no damping in the system. We note that Pauli block- 
ing may be suppressed when the double exciting state is 
formed by electron/hole pairs with opposite spins. 

Since the nonlinear signal vanishes for ideal bosons, 
one can recast it to the alternative simplified form as if 
from the ESA from otherwise Pauli blocked N states as 



N I,, .|2|,, .|2 



(50) 



■ ^ (-wi - E*) (lj 3 - Ei ~ E j + E$) 
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At zero time delay T2 = we have only the diagonal res- 
onances. The states with small damping dominate the 
picture. However as the time delay progresses the off 
diagonal resonances appear, as demonstrated in Fig. [6] 
It is convenient to interpret the signal by comparing it 
with the linear absorption (the top-marginal graph in the 
figure). For our numerical simulations, we chose the po- 
tential kink height to be V\ — —V2 — 20. Unless stated 
otherwise, all the energies are in units of 3joa/R. We 
have also added constant dephasing rate 7 = 0.1 to ac- 
count for possible contact with an external phonon bath. 
We have also limited ourselves to electronic states with 
angular momentum up to m = 9/2. An idealized single 
graphene layer grown epitaxially on SiO was considered. 

The tallest absorption peak comes from the transition 
5/2 — !> 7/2, and is located at u>i rj 40.0. That signi- 
fies the true bound state when the electron (hole) energy 
reaches the height of the potential barrier (See Eq. ( 29 ) 
and Fig. |3]). The remaining absorption peaks correspond 
to quasi-bound states with finite lifetime. For the linear 
spectrum, the latter brings the peak broadening. To ex- 
tract additional information about the dynamics of the 
quasi-bound states, we resort to the photon echo signal. 
At zero time delay t-i — 0, this provides the same in- 
formation as the linear absorption. The positions and 
magnitudes of the cross-peaks (rapid change in the sign 
of the signal) on the main diagonal correspond to those in 
the linear absorption. The existence of the signal comes 
from the fact that the electrons (holes) are not coupled to 
a simple bath of harmonic oscillators (constant dephas- 
ing) . The pattern of the cross-resonances along the main 
diagonal is the manifestation of the destructive interfer- 
ence between the GSB, ESE on one side and the ESA 
pathways on the other. The latter takes into account 
the Pauli blocking effect on the biexciton (two electron- 
hole pair) states. At this point, we completely neglected 
the effect due to the Coulomb interaction between elec- 
trons. Later, we shall demonstrate that it is a reasonable 
assumption for small QDs. Thanks to the very simple 
exciton scattering matrix based on Pauli blocking, only 
(Eq. (44)), we can employ a simplified quasi-particle pic- 



ture in order to describe the signal (Eq. (50)) 



By increasing the time delay T2, we may monitor the 
lifetime dynamics of the quasi-bound states as follows. 
The ESA and ESE contributions to the signal is reduced 
and finally the GSB signal survives (the lowest graphs 
in Fig. [6]). In between, the off-diagonal cross-peaks ap- 
pear at a time. Those with the smaller dephasing rate 
(strongly bound to the QD) appear first. The most pro- 
nounced cross-peaks are those which are correlated to the 
true bound state. 

We next turn our attention to the dressed Dirac elec- 
trons confined to the potential induced QD. The dressing 
opens up a dynamical gap which can be controlled by the 
intensity and polarization degree of CW pumping light. 
We shall probe the dynamic gap by the photon echo tech- 
nique described above and compare it to the linear ab- 
sorption. The gap allows for many more bound states 
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FIG. 6. Photon echo signal for various time delays T2 calcu- 
lated from Eq. (46). The left panels correspond to E g — 0, 



while the right panels demonstrate the dressing effect E g 
10. The top-marginal graphs show the linear absorption. 
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since the wave vector of the outgoing electronic wave k 
can cross over into the imaginary plane, thus effectively 
quenching the outgoing wave and bounding the electron 
states (See Fig. [3]). For our simulations, we chose the gap 
Eg = 10 which may be achieved either for small QD or 
intense pumping field with circular polarization. We note 
that the gap may also be induced by a polar substrate^. 
The gap achieves several bound states for the 1/2 — s- 3/2 
electronic transitions. Since the wave functions for larger 
angular momentum are highly oscillatory, the latter tran- 
sition posses highest oscillator strength, thereby effec- 
tively shifting the position of the main peak in the linear 
absorption (see Fig. [6] right panel) . The remainder of 
the peaks also contain a mixture between the bound and 
quasi-bound states. To separate these, we shall look at 
the photon echo at r 2 > 0. Finally, the resulting GSB 
reveals the truly bound states (see Fig. fright panel). 

To examine the role played by Coulomb scattering, 
we shall employ the full form of the scattering matrix. 
The approach is based on the nonlinear exciton equations 
(NEE). We refer the reader to the comprehensive review 
of the technique given by Abramavicus and Mukame^. 
Exciton scattering is best described in the eigenstate ba- 
sis of Eqs. (28) and (27). Keeping in mind that we 



can have at most two excitons, leads to effective trunca- 
tion schematics of otherwise infinite series of intertwined 
NEEs 1 ^. In the latter case, an appropriate factorization 
scheme has been applied. We have also neglected inco- 
herent exciton transport. 

The photon echo and the linear absorption are shown 
in the right panel of Fig. [7] for larger size QD. We see the 
off-diagonal correlation resonances and symmetry break- 
ing for T2 = 0. These indicates the bonding and anti- 
bonding biexciton resonances with the biexciton binding 
energy of a few eV. Indeed, when the biexciton binding 
energy is increased as a result of the Coulomb interac- 
tion, the ESA peaks are shifted along lu^: downwards for 
positive anti-binding (exciton repulsion) and upwards for 
negative bonding energy (exciton attraction). The ESA 
cross-peaks are no longer cancelled by the GSB and ESE, 
thus creating the doublets. By increasing the QD size, we 
see the formation of excitons with exciton binding energy 
of 10 — 30 eV in the left panel of Fig. [7J Signatures of 
the off-diagonal quadratic coupling also persist for longer 
delay times r 2 > 0. 



VII. CONCLUDING REMARKS 

We proposed dressing the Dirac electrons with circu- 
larly polarized photons in order to localize them within a 
QD on graphene monolayer. We also investigated the lo- 
calization of dressed electrons in a cylindrical QD formed 
on bilayer graphene. When graphene is irradiated with 
a circularly polarized electromagnetic field, An energy 
gap opens up in the dispersion relation for graphene in 
the presence of this electromagnetic field. Consequently, 
the resulting confined electronic states for a QD seem 



3, 




3 °" 3 

~3 0.2 

So 0.1 

°° 0.0 



34 36 38 40 42 44 46 ' 

Sf 1 '(a>3»0.,-a>i) 



34 36 38 40 42 44 46 . 
SfV 3 ,0. ,-CJi) 




34 36 38 40 42 44 46 48 




34 36 38 40 42 44 46 . 



Sf(a) 3 ,0.2,-o>{) 



48 






46 




j& - 


44 






«, 42 
3 40 










38 
36 






34 







34 36 38 40 42 44 46 48 
-0) l 

Sf ) (w 3 ,0.8,-wi) 



34 36 38 40 42 44 46 48 



syV 3> o.8,-wi) 




Ml 













34 36 38 40 42 44 46 48 
-£<Jl 

$"(.0)3,1.5 ,-COi) 



34 36 38 40 42 44 46 48 




34 36 38 40 42 44 46 48 
-0)i 



34 36 38 40 42 44 46 48 



FIG. 7. The same as in Fig. [6] w ith the Coulomb scattering 
taken into account via Eq. ( |45| with the scattering matrix 
given by Eq. The left panels represent large QD, R/ao = 
1000. Right panels are for smaller QD, R/ao = 100. 
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to have properties that are similar in nature to the sur- 
face states of topological insulators. Their energies are 
located inside the energy gap and the wave functions de- 
cay as a function of distance from the interface of the 
potential. These topological states are robust with re- 
spect to the effects of disorder. Our calculations showed 
that the dressing does not only open a dynamical gap 
in the energy dispersion spectrum, but it also leads to a 
renormalization of the Fermi velocity as well as the in- 
tra layer and inter layer coupling parameters. In fact, 
in the bilayer configuration, the dressing serves as a tool 
for tuning the energy gap. That is, it can either close 
or open the gap, depending on the polarity of the po- 
tential and the direction of the light polarization. Linear 
spectroscopy cannot resolve the short lived broadened 
excitonic states and must be resolved by using a four- 
wave mixing technique known as photon-echo. This elim- 
inates the inhomogeneous broadening due to impurities, 
and to focuses on the intrinsic lifetimes of the electronic 
states. We measure the localization through the elec- 
tronic density of states, The strong dynamical screening 
of the Coulomb interaction leads us to consider only the 
Pauli blocking due to the Fermi statistics. We simplify 
the signal interpretation by switching to the quasipar- 
ticle picture. Those are give as the deviation from the 
harmonic oscillator for which the nonlinear signals dis- 
appear. This allows us to consider only excited states 
absorption Liouville pathways. In this way, we are able 
to reduce the interference due to the usual combination 
between the ground state bleaching and excited states 
emission. Visible light is used to map the QD interband 
transitions onto 2D spectra and terahertz pulse shaped 
fields for the intraband transitions. Important aspects 
of terahertz pulse shaped fields for the intraband transi- 
tions will be reported elsewhere. The latter will allow us 
to use a novel and more convenient phase cycling method 
to obtain the response 1 ^. 



ACKNOWLEDGMENTS 



4 G. Pal, W. Apel, and L. Schweitzer, "Electric transport through 

circular graphene quantum dots: presence of disorder," Arxiv 

preprint arXiv:1107.0113(2011). 
5 C. H. Park and S. G. Louie, "Tunable excitons in biased bilayer 

graphene," Nano letters 10, 426-431 (2010). 
6 A. Matulis and F. M. Peeters, "Quasibound states of quantum 

dots in single and bilayer graphene," Physical Review B 77, 

115423 (2008). 

7 B. Wunsch, T. Stauber, and F. Guinea, "Electron-electron inter- 
actions and charging effects in graphene quantum dots," Physical 
Review B 77, 035316 (2008). 

8 P. Hewageegana and V. Apalkov, "Electron localization in 
graphene quantum dots," Physical Review B 77, 245426 (2008). 

9 0. V. Kibis, "Metal-insulator transition in graphene induced 
by circularly polarized photons," Physical Review B 81, 165433 
(2010). 

10 Y. Zhou and M. W. Wu, "Optical response of graphene under 
intense terahertz fields," Physical Review B 83, 245436 (2011). 

11 D. Abramavicius, B. Palmieri, D.V. Voronine, F. Sanda, and 
S. Mukamel, "Coherent multidimensional optical spectroscopy 
of excitons in molecular aggregates; quasiparticle versus super- 
molecule perspectives," Chemical reviews 109, 2350-2408 (2009). 

12 W. Kuehn, K. Reimann, M. Woerner, T. Elsaesser, R. Hey, and 
U. Schade, "Strong correlation of electronic and lattice excita- 
tions in gaas/algaas semiconductor quantum wells revealed by 
two-dimensional terahertz spectroscopy," Physical Review Let- 
ters 107, 67401 (2011). 

13 E. G. Kavousanaki, O. Roslyak, and S. Mukamel, "Probing ex- 
citons and biexcitons in coupled quantum dots by coherent two- 
dimensional optical spectroscopy," Physical Review B 79, 155324 
(2009). 

14 C. Gerry and P. Knight, Introductory Quantum Optics (Cam- 
bridge Univercity Press, 2003). 

15 S. Mukamel, Principles of Nonlinear Optics and Spectroscopy 
(Oxford University Press, 1995). 

16 S. Mukamel, R. Oszwaldowski, and D. Abramavicius, "Sum-over- 
states versus quasiparticle pictures of coherent correlation spec- 
troscopy of excitons in semiconductors: Femtosecond analogs of 
multidimensional nmr," Physical Review B 75, 245305 (2007). 

17 V. Chernyak, W. M. Zhang, and S. Mukamel, "Multidimensional 
femtosecond spectroscopies of molecular aggregates and semicon- 
ductor nanostructures: The nonlinear exciton equations," The 
Journal of chemical physics 109, 9587 (1998). 

18 The interaction is modified by the dressing, as indicated by tilde. 

19 V. M. Agranovich, "Electrodynamics of excitons in two- 
dimensional systems," (Taylor & Francis, 1992). 

20 A. Qaiumzadch and R. Asgari, "Ground-state properties of 
gapped graphene using the random phase approximation," Phys- 
ical Review B 79, 075414 (2009). 



The authors gratefully acknowledge the support of 
Air Force Research Lab (AFRL) by contract No.# FA 
9453-11-01-0263; the National Science Foundation (NSF) 
through Grant No.# CHE-1058791, DARPA BAA-10-40 
QuBE; from Chemical Sciences, Geosciences, and Bio- 
sciences Division, Office of Basic Energy Sciences, Office 
of Science, (U.S.) Department of Energy (DOE). 



Appendix A: Derivation of Eq. |T|) 



1 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, 
and A. K. Geim, "The electronic properties of graphene," Re- 
views of modern physics 81, 109-162 (2009). 

2 L. J. P. Xavier, J. M. Pereira, A. Chaves, G. A. Farias, and F. M. 
Peeters, "Topological confinement in graphene bilayer quantum 
rings," Applied Physics Letters 96, 212108 (2010). 

3 T. Paananen, R. Egger, and H. Siedentop, "Signatures of wigncr 
molecule formation in interacting dirac fermion quantum dots," 
Physical Review B 83, 085409 (2011). 



The original Hamiltonian has the form 



n = Hv F (T- (k- —A 
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The vector potential operator of the electromagnetic field 
can be partitioned as: 



A = A + A 

i=l 

l2irhc 2 ( 



e + ao + e_a ( 



(A2) 
(A3) 
(A4) 



Here y2e± = e x ± ie y are polarization vectors given 
in terms of the unit vectors along corresponding Carte- 
sian directions; £1 is the mode quantization volume. As 



one can see Eq. ( A2 ) describes the electromagnetic wave 
propagating along z— axis (transverse to graphene). It 
is clock-wise circularly polarized. We will need the 
circular polarization since graphene is gapless and no 
RWA is applicable. The rest of the optical modes de- 
scribed by Eq. (A3) are linearly polarized. Note that 



we have no phase on the optical filed since we assume 
graphene being ideally flat and situated at z = 0. That 
is e-xp(±ik z z) = 1. Substituting Eq. ( |A2"| ) into Eq. (KJ\ 
and denoting — e^inhvp /ujifl = Wi/y/Nj, in order to 
keep notation consistent with RefP^ we obtain Eq. ([!]) . 



Appendix B: Derivation of Eqs. (10) and (111 



We first need the following identities 



ff±|±,JV ) = |T,iV > 
<t±\T,N ) =0 

Therefore, we shall have: 



(Bl) 
(B2) 



Hv F (a x k x + ayk y ) |V>±,iv ) (B3) 
= hv F ((cr + + oJ)k x + i(cr_ - (T+)k y ) \ip±, Na ) 
= hv F + (T-)k x + i{(J- — <J+)k y ) 
x (cos </>| ±, N ) ± sin 0|=p, JV ± 1)) 
= hv F {k x (cos0|T,iVo)±sm0|±,iVo±l))± 
±ik y (cos0|=F, A ) T sin0|±, N ± 1))) 



Using the above equation we can calculate all the nec- 



essary matrix elements: 

H x =H l + V{x,y) 

{^+,N \'fLl\^+ l N ) 
= hv F ((cos(j)(+,No\ + sin (/)(-, A + 1|) 
xk x (cos0|-, N ) + sin0|+,7V o + 1)) 
+ik y (cos0|-, A ) - sin0|+, N + 1))) = 

(l/->-,N \Til\lP-M ) 

= Hv F ((cos<f>(-,N \ - sin0(+, A - 1|) 
xk x (cos0|+, A ) - sin0|-,7V o - 1)) 
-iky (cos0|+, N ) + sin0|-, N - 1))) = 

= hv F ((cos 0(-, Nq\ - sin </)(+, A - 1|) 
xk x (cos 0|- N ) + sin0|+, N + 1)) 
+ik y (cos0|-,7V o ) - sin0|+, N + 1))) 
= cos 2 (k x + iky) 

(<0+,JV o |^x|'0-,JVo) 

= hv F ((cos(j)(+,No\ +sin0(-, N a + 1|) 
xk x (cos0|+,A o )-sin0|-,7V o -l>) 
-iky (cos0|+,7V o ) + sin0|-,A o - 1))) 
= cos 2 (fe x — iky) 
(ip±.N \V(x,y)\ip ±tNo ) 
= V(x, y) (cos 2 + sin 2 0) = V(x, y) 



(B4) 



(B5) 



(B6) 



(B7) 



(B8) 



For 'H.-z matrix elements we will need, the following 
identities. 

(^±,A r ol (J + +O'-|'0±,JVo) (B9) 

= (*0±,iVo I (cos 0|T, A ) ± sin 0|±, N ± 1) ) 
- (cos0(±,7vo|±sin0(T,A o ±l|) 
x (cos0|=F, No) ±sin0|±, A ± 1>) = 

(ip TtNo \a + + a-\tp± tNo ) = cos 2 (BIO) 



Appendix C: Derivation of Eqs. (19 1 and (20) 



For the b-layer we will need the following identities: 

(■^N \<J±\lpN ) (CI) 

_ ( (^+MM±\^+M ) (^+,N W±\^-,N„) 
(^-,N W±\i>+,N ) (lp-,N \<r±\4>-,No) 
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FIG. 8. Feynman diagrams for the double quantum coher- 
ence signal generated in k s = ki + k2 — k3 direction. 



Appendix D: Double quantum coherence. 

The double quantum coherence signal can be derived 
from the diagrams in Fig. [8] assuming the following form: 

^f^ii*(«3,^.n = Q) (Di) 

= 2lm ^ Mel'Ve2 j2 Me3^4 G 'e4(w3)Gg3(w2 - U) 3 ) 
el,e2,e3,e4 

x [T'e4,e3;e2,el(k'3 + E e3 )G e 2^el (^3 + E e {) — 
— r e 4 : e3;e2,el (w2)G e 2,el (^2)] 

When the Coulomb scattering may be neglected the 
above signal is greatly simplified into sum-over-states ex- 
pression with the explicit Pauli blocking principle: 



CJ' 1 J'2j"3,j4 / \ 
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This signal vanishes identically despite the Pauli induced 
scattering making it a measure of the screened Coulomb 
interaction. 



